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ANGLE-OF -ATTACK FLOW ABOUT AXISYMMETRIC 
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By Richard W. Barnwell 
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SUMMARY 

A time-dependent numerical method for calculating supersonic flow about blunt 
bodies at large angles of attack is presented. The axisymmetric bodies with sharp shoul- 
ders which are treated are constructed with a generator composed of segments of constant 
curvature. The nonaxisymmetric bodies have continuous slope and curvature. All flow 
fields are inviscid and adiabatic and have one plane of symmetry. 

A modification to the method of characteristics is introduced for use at the shock 
wave. A two-step finite-difference method of second-order accuracy is used at the body 
surface and in the region between the shock and body. A new finite-difference technique 
is introduced for use at sharp sonic shoulders. 

Comparisons of the results of the present method with experiment and the results of 
other methods are made for the flow of equilibrium air past the Apollo command module 
at the trim angle of attack and for perfect gas flow past a spherical cap and a spherically 
blunted cone at angle of attack. Both the cap and the blunted cone are terminated with 
sharp shoulders. Results are presented also for perfect gas flow past a prolate spheroid 
with its major axis normal to the flow. 

INTRODUCTION 

Time-dependent finite-difference methods provide a means of treating the problem 
of inviscid supersonic flow past a blunt body as an initial-value problem since the equa- 
tions for inviscid transient flow are always hyperbolic. Results for steady flow are 
obtained from the asymptotic solution to the transient problem. One of the major advan- 
tages of these methods is that there is no conceptual difficulty in extending them to treat 
such three-dimensional effects as angle of attack and nonaxisymmetric body geometry. 

In general, the chief difficulty encountered in making a three-dimensional, rather than a 



two-dimensional, time-dependent calculation is the additional time required to perform 
the computation. 

Time-dependent methods for calculating three-dimensional blunt-body flow fields 
have been developed by Rusanov (ref. 1), Bohachevsky and Mates (ref. 2), Moretti and 
Bleich (ref. 3), and Xerikos and Anderson (ref. 4). Only the method of reference 1 can be 
applied to anything but axisymmetric bodies. All the methods except that of Moretti and 
Bleich produce results of first-order accuracy in the mesh spacings; the method of ref- 
erence 3 produces answers of second-order accuracy. The method of Bohachevsky and 
Mates requires a much larger number of grid points than the other methods because the 
bow shock wave is treated as an internal feature of the flow rather than as a boundary. 

As a result, a great deal more computer time is required for this method than for the 
others. 

Cohen, Foster, and Dowty (ref. 5) have employed the refined Godunov method of 
Masson, Taylor, and Foster (ref. 6) to develop an approximate time-dependent method for 
calculating angle-of-attack flow. The flow is calculated in the plane of symmetry, and 
trigonometric functions are used to approximate the cross-flow derivatives. 

The purpose of this paper is to present a time -dependent method for calculating 
three-dimensional flow fields about two fairly general classes of bodies. The first class 
is that of axisymmetric bodies and includes bodies with discontinuous surface slope and 
curvature. The second class is that of nonaxisymmetric bodies with one plane of sym- 
metry and with continuous surface slope and curvature. Both perfect and equilibrium gas 
models can be treated. A previous version of this method for calculating flow about axi- 
symmetric bodies with sharp shoulders at angle of attack was presented in reference 7. 
The present method is a refinement of a previous method presented in references 8 and 9 
and extends the applicability of that method to three dimensions. As was the case for the 
method of references 8 and 9, the present method is closely related to that of Moretti and 
coauthors (refs. 3 and 10). 


SYMBOLS 


A,B,C,CO 
D,E,E',F J 


points in figures 4 and 5 


a,b,c 


matrices defined by equations (B2) 


Ai,Bi,Ci7) 
Di,Ei J 


quantities defined by equation (6) 
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speed of sound 


semimajor and semiminor axes of prolate spheroid 

quantities defined by equations (15) 

vector defined by equations (B2) 

matrix defined by equation (Bll) 

internal energy 

eigenvalue of matrix E 

quantity defined by equation (B17) 

unit vector normal to shock, defined by equation (26) 

unit vector in free-stream direction, defined by equation (31) 

unit vectors in x-, y-, and (^-directions, respectively 

quantities used in equation (18) 

lower and upper bounds for inequality (B20) 

quantity defined by equation (B24) 

quantity defined by equation (B13) 

maximum value of F u (x,e) for given value of e 

quantity defined by equation (25) 

matrix defined by equation (B12) 

eigenvalue of matrix G 


total enthalpy 



h 

h y >h ?J 

K 

Lx>Ly ,L^J 

M 

m 

m 


P1.P3 

P 

Po 

Pt 

Q 

R 

r 

r n 

r l’ r 2’ r 3 

s 

t 

U,V 


static enthalpy 

direction cosines of bicharacteristic in t,x,y,<p space 
curvature of reference surface 

wavelength of error solutions in x-, y-, and (^-directions, respectively 
Mach number 

quantity defined by equations (B15) 
quantity defined by equations (B22) 
quantities defined by equations (B22) 
pressure 

reference pressure, 1 atmosphere (101.3 kN/m2) 

exact value of stagnation pressure 

quantity defined by equation (B23) 

quantity equal to right side of equation (33) or (35) 

perpendicular distance from coordinate axis 

nose radius 

quantities defined by equations (B8) 

distance along surface from axis in plane of constant cp 
time 

components of velocity tangent and normal to shock, respectively 
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Ux>Uy,U,^ 

u 


u,v,w 

u,v 

V 

V 

W 

X,x 

V 


y 


components of U in x-, y-, and (^-directions, respectively 
quantity defined by equations (B15) 

components of velocity in x-, y-, and (^-directions, respectively 

velocity components shown in figure 5 

magnitude of total velocity 

quantity defined by equation (B3) 

vector defined by equations (B2) 

distance along generator of reference surface from axis 
normalized coordinate defined by equations (7) 
coordinate along normals to reference surface 


y = y - Yb 


z 

a 

a 


Q k’ Q y > a <p 


y 

y 


coordinate along axis 
angle of attack 

exponent used in equation (B6) 
quantities defined by equations (30) 
quantities defined by equations (24) 
ratio of specific heats 

ratio of static enthalpy to internal energy, h/e 


at,ax7i 
AY,A$ J 


mesh spacings for time r and coordinates X, Y, and 


5 


shock layer thickness, function of t,X,# 


6 

e damping coefficient 

e = — - — 
e max 


C 




e 

x 


v 


v vV<2> 


P 


a 


7 


$,<p 


X 

w 

Subscripts: 


distance from corner, x c m ^ n - x 

quantities defined by equations (B7) 

angle between normal to reference surface and axis 

scale factor for x-coordinate defined by equation (3) 

quantities defined by equations (B15) 

exponent employed in equation (18) 

quantities defined by equations (B22) 

density 

angle between normal to shock and free-stream direction 
time 

azimuthal angle 

quantity which varies from 0 to 1 
small nonnegative number 


b properties at body 

c properties at shoulder 

I initial solution 
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max 


maximum value 


min minimum value 

s properties at shock 

°° properties in free stream 

90° properties for cp = 90° 

Superscripts : 

o> exponent between 0 and 1 

* conditions when u = a and v = 0 

• differentiation with respect to t or t of functions which do not depend 

on Y 

' differentiation with respect to x or X of functions which do not depend 

on Y 

differentiation with respect to cp or <i> of functions which do not depend 
on Y 

— vector quantity 


ANALYSIS 

The present method for calculating numerical solutions for time-dependent, inviscid, 
three-dimensional flow past blunt bodies traveling at supersonic speeds is described in 
this section. A number of grid points are located on the bow shock wave, the body sur- 
face, and between the shock and surface. The region of computation must contain the 
entire zone of subsonic flow. An initial solution, which can be quite general, is assumed; 
and the flow at each of the grid points is calculated for a number of time steps. At each 
cycle of the computation an initial-value problem is solved to determine the solution at 
the new time step from the solution at the previous time step, subject to the appropriate 
boundary conditions. Results for steady flow are obtained after many time steps when 
the time derivatives of the flow properties are sufficiently small. It should be noted that 
the locations of the grid points adjust with time as the location of the bow shock adjusts. 
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A two-step finite-difference approximation to the time-dependent method of character- 
istics is used at the bow shock wave, whereas the two-step, time-dependent, finite- 
difference scheme of Brailovskaya (ref. 11) is used at the surface and between the shock 
and surface. The solution for a given time step is determined first at the points on the 
shock wave, then at those between the shock and body, and finally at those on the body 
surface. 

Both the present method and the method of references 8 and 9 can be used to calcu- 
late axisymmetric flow about blunt bodies with sharp sonic shoulders. A major advantage 
of the present method is that it is not necessary to specify any of the flow properties at 
the shoulder. With the method of references 8 and 9, it is necessary to specify that the 
Mach number at the shoulder is 1 when the flow upstream of the shoulder is subsonic. A 
second advantage of the present method is that the computational techniques which are 
used at the bow shock wave and body surface are much more efficient than the conven- 
tional time-dependent method of characteristics which is used in references 8 and 9. 


Basic Coordinate System and Governing Equations 

The basic coordinate system which is used is similar to that employed in refer- 
ences 8 and 9. An axisymmetric reference surface is established as shown in figure 1. 
The coordinates are the azimuthal angle cp and the distances x and y, which are 
measured along the reference surface and normal to it, respectively, in planes of con- 
stant cp. The components of velocity in the x-, y-, and ^-directions are u, v, and w. 
The angle in a plane of constant cp between the normal to the reference surface and the 
direction of the coordinate axis is designated as 9 and satisfies the differential equation 



( 1 ) 


where K is the local curvature of the reference surface. The distance r from the 
axis is given by the equation 


pX=X 

r(x,y) = \ cos 0(x)dx + y sin 0(x) 
J x=0 


As stated previously, the bodies which the present method will treat have one plane 
of symmetry. The free-stream velocity vector is parallel to this plane. The angle <p 
is measured from the leeward side of this plane, and the reference surface is constructed 
of segments of constant curvature. The reference surface shown in figure 1 is constructed 
of one segment and hence is spherical. 

The type of coordinate system to be used with an axisymmetric body with a sharp 
shoulder is shown in figure 2. The reference surface which is used with the body in the 
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figure has three segments of constant curvature. The coordinate system is focused at 
the corner in the manner of references 8 and 9 so that in the vicinity of the shoulder 
(^ x c,min = x = x c,max) line °* intersection of the reference surface and a plane of con- 
stant cp is a circular arc centered at the shoulder. The use of a coordinate system 
which focuses at sharp shoulders is crucial to the present method. 


The coordinate system to be used with an axisymmetric body without sharp shoul- 
ders but with a generator composed of segments of constant curvature is similar to that 
used with bodies with sharp shoulders in that the reference surface is located at a con- 
stant distance y^ from the body along the normals to these surfaces. For nonaxisym- 
metric bodies, the distance y^ is generally a function of x and cp. 

The equations for the conservation of mass and momentum, written in dimensional 
form in terms of the time t and the coordinates x, y, and cp, are 


Continuity: 


dp 

dt 


+ P 


1 8u 
X 0X 



1 9W 
r Bcp 


cos 0 (k sin 

+ u + — + v 

r \X r / 


= 0 


(2a) 


x -momentum: 


— + — -^2 + — uv - £25LJl w 2 _ o (2b) 

dt Xp 9x X r 


y-momentum: 


dv 1 8 P K 9 

dt P 0 y X 


sin 9 o n 
w z = 0 


(^-momentum: 


(2c) 


dw 

dt 


1 0 p cos 9 


uw + 


sin 9 
r 


vw = 


0 


(2d) 


where p and p are the pressure and density of the gas, X is the scale factor for the 
x-coordinate and satisfies the equation 


X = 1 + yK (3) 

and d/dt is the total derivative with respect to time along a streamline and can be 
written as 
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d _ 9 [ u 9 , y 8 W 9 
dt 8t X 9 x + 9 y + r 8 <p 


In this paper, the energy equation is used in the forms 
de _ _p_ dp = o 


dt p2 dt 


( 4 a) 


and 


^P + a 2 p 

dt 


_1 jHa + 9v + j._9w + cos 0 u + (K + sin 0 ' 
X 9x 9y r 8<p r \x r 


= 0 


( 4 b) 


where e and a are the internal energy and speed of sound of the gas, respectively. 

It should be noted that equations (2) and (4a) can be written in conservation form as 


9Ai 9Bi 9Ci 9Di 

— i + — - + — i + — i + E i = 0 

9t 9x 8y 9 cp 


where 


Ej = (pu cos 0 + pv sin 0)£ 

E2 = p( u 2 - w2)A cos q + p U v^K + p sin 0j 

E3 = puv cos 0 + p(v 2 - w 2 )^ sin 0 - K(p + pu 2 ) 

E4 = 2(puw cos 0 + pv w sin 0)^ 


E5 = (puH cos 0 + pvH sin 0)£ 


(i = 1,2, 3, 4, 5) 


Ai = PA. 

Bi = pu 

Cj = pvX 

Dj = pw y 

A2 = puX 

B2 = p + pu 2 

C2 = puvX 

D2 = puw j 

> 

CO 
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B3 = puv 

C3 = (p + pv 2 )X 

D3 pv w£ 

A4 = pwX 

B4 = puw 

C4 = pv wX 

D4 = (p + pv 

A5 = (pH - p)X 

B5 = puH 

C5 = pvHX 

D5 = pwH y 


■x 


( 5 ) 


( 6 ) 


y 


10 



The quantity H in equations (6) is the total enthalpy and is written as 


H = h + -i(u2 + v2 + w2) 

£t 

where h is the static enthalpy of the gas. 

Computational Coordinates 

The shock-wave and body-surface locations are specified by the equations 


y = y s (t,x,<p) 


and 


y = yb( x ><?) 


respectively, and the distance between the shock and surface for given values of t, x, 
and cp is 6 = y s - y^. The function y^ for the body surface is known, and the func- 
tions y s and 6 must be determined as part of the solution. 

In this paper, calculations are made in terms of the independent variables 


r = t X = x 


Y = 


y-yb 

6 


= cp 


( 7 ) 


It should be noted that the geometrical variables X, Y, and $ do not form an orthog- 
onal set. The partial derivatives with respect to t, x, y, and cp are related to those 
with respect to t, X, Y, and <3? as follows: 


8 _ 8 Y6 9 
8t ~ 8 t " 6 8Y 


8 

8x 

8 

9y 


— - -(y 6' + yfA— 

8X S\ ^ b /8Y 

1 8 

6 8Y 


(8) 


— — = — — - Y6 + y h — 

8 cp 8$ S'' D /8Y 

J 

where 6, 6*, and 6 are the derivatives of 6 with respect to r, X, and 

respectively. 
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The shock layer to be calculated is subdivided with the X-, Y-, and <&- coordinate 
surfaces as shown in figure 3. The planes of constant <£ are spaced from $ = 0° to 
$ = 180° with a uniform mesh spacing A<$. Since the flow fields which are treated in 
this paper have one plane of symmetry, only one side needs to be calculated. The coor- 
dinate Y has values of 0 and 1 at the body and shock, respectively. A uniform mesh 
spacing AY separates the surfaces of constant Y which are located between the sur- 
face and shock wave. The coordinate X is measured along the generator of the axisym- 
metric reference surface. As stated previously, this generator is constructed of seg- 
ments of constant curvature (the generator illustrated in fig. 3 has two such segments). 
The surfaces of constant X are spaced so that the distance separating them along the 
generator of the reference surface is a constant spacing AX. It should be noted that the 
surfaces of constant X are orthogonal to the reference surface. 

The flow properties are calculated at the intersections of the surfaces of constant 
X, Y, and <£. Since the Y-coordinate is normalized with respect to the distance between 
the shock wave and the body surface, the locations of the grid intersections move as the 
shock wave adjusts to its steady position. 

Calculation of Flow Within Shock Layer 

The procedure for the calculation of flow within the shock layer is used at the grid 
points between the shock and surface. A related approach is used at most of the points 
on the surface. 

General procedure .- The equations which are solved are equations (5) in terms of 
the independent variables r, X, Y, and $ and are written as 


(i = 1,2, 3, 4, 5) (9) 

The quantities y^ and y^ are the derivatives of y b with respect to X and 4>, 
respectively. 

The equations are used in this particular conservation form because the conserved 
functions have continuous derivatives in the vicinity of sharp sonic shoulders although 
some of the flow properties have unbounded derivatives at these shoulders. The advan- 
tages of using this form to obtain solutions on streamlines near sharp shoulders were 
discussed in reference 9, and the advantages of using it on surface streamlines are dis- 
cussed subsequently. 
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9X 6 9Y 


Ci - Y5Ai - (y(j + Yfi')Bi - (y b + y6)d 


^H- El+ i(6A 1 + 6 'B i + SD i )J> 
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As stated in reference 9, the derivatives with respect to Y in equations (9) are 
well behaved at points in the flow field near sharp shoulders because the quantities which 
are differentiated are proportional to Y at these shoulders. This follows because the 
segment of the generator of the reference surface associated with the shoulder is a cir- 
cular arc with its origin at the shoulder as shown in figure 2. Hence, is constant at 
the shoulder so that 


y{> - 0 Yb - 0 -^c,min = x = ^c,max 

The scale factor X at points near the shoulder satisfies the equation 



Note that the quantity y^ is always negative for a sharp shoulder. Since the quantities 
q are all proportional to X, it follows that these terms are proportional to Y near 
the shoulder. In view of these considerations, the quantities which are differentiated 
with respect to Y in equations (9) satisfy the following relationship near sharp 
shoulders: 


Ci - Y6Ai - (Y5 1 + y^Bi - (YS + y^D^ = q - Y(6Ai + 6’Bi + 6Di) °c Y 

A two-step, time-dependent, finite-difference method of second-order accuracy in 
the mesh spacings AX, AY, and A# and first-order accuracy in At is used to inte- 
grate the governing equations. This method was first used by Brailovskaya (ref. 11). In 
the first step the time derivative is approximated with a forward-difference expression, 
whereas in the second step the time derivative is replaced with a backward-difference 
expression. The first-step solution, which is designated with the subscript I, is deter- 
mined with the equations 

A 1>i (t,X,Y,*) = Ai(T-Ar,X,Y,$) + At ±- Ai(T-AT,X,Y,4>) (i = 1,2, 3, 4, 5) (10) 


The partial derivatives SAi/W in equations (10) are determined by evaluating the right 
side of equations (9) at the point (t - At,X,Y, 4>). Unless otherwise noted, the partial 
derivatives with respect to X, Y, and 4> in equations (9) are approximated with 
central-difference formulas. For example, the partial derivatives 9Bt/9X in equa- 
tions (9) are approximated with the expressions 
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— Bi(T-AT,X,Y,$) = — ir-rBi(T-Ar,X+AX,Y,$) - B i (r-Ar,X-AX,Y,$) 
0X 2 AX _ 


The second-step solution is obtained with the equations 


Ai(r,X,Y,*) * A i (T-Ar,X,Y,#) + At ^ A i;I (T,X,Y,$) - e <^(AX) 4 Ai(r-AT,X,Y,*) 
+ (AY) 4 Ai(T-AT,X,Y,$) + (AS) 4 — Ai(T-AT,X,Y,$)} 

0Y 4 a# 4 J 


(i = 1,2, 3, 4, 5) (11) 

where the values of the partial derivatives SA^j/st at the point (r,X,Y,$) are deter- 
mined by evaluating the right side of equations (9) at this point with the first-step solu- 
tion. The terms of fourth order in equations (11) are nonphysical damping functions 
similar to that used by Richtmyer and Morton (ref. 12) which are added to eliminate 
instabilities. Values for these derivatives are determined with five-point formulas of 
the form 

, s 4 

(AX) 4 — Ai(T-AT,X,Y,$) = A i (r-AT,X+2AX,Y,4>) + Ai(r-Ar,X-2AX,Y,4>) 

8X 4 

- 4jAi(T-AT,X+AX,Y,<f») + Ai(T-Ar,X-AX,Y,4>) 

+ 6Ai(T-AT,X,Y,4>) 


Five-point formulas of this type cannot be used to evaluate the terms (AY) 4 


9 4 Ai 


9Y 


on 


surfaces of constant Y adjacent to the shock wave and body surface because points with 
the coordinates Y + 2 AY would be beyond the shock and those with coordinates 
Y - 2 AY would be inside the body. On the surfaces of constant Y next to the shock 


wave and body surface, the terms (AY)^ 


a 4 Ai 

0Y 4 


are replaced with terms of the form 
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8 3 Ai . . 9 9 3 A 


-(AY) 3 and (AY) 3 -, respectively. These quantities are evaluated with the 


0Y*- 


8Y‘ 


following four-point formulas: 


8^ A* 

±(AY) 3 = A i (r-AT,X,Y±2AY,<i>) - Ai(T-AT,X,Y*AY,$) 

9Y 3 


- 3 Ai(r-AT,X,Y±AY,$) - Ai(T-AT,X,Y,‘f>) 


The density and the velocity components are determined from the quantities Aj, 
A2, A3, and A4 with the equations 


P 





w 


A4 

Al 


The value of the scale factor X at time t, which is needed in order to determine the 
density, is calculated with the following equation obtained by substituting the third of 
equations (7) into equation (3): 


X = 1 + (6Y + y b )K 


Since it is necessary to know the shock layer thickness 6 at time t, the solution at the 
shock is determined first at each time step. The internal energy of the gas is determined 
with the equation 


e = — - - — (u^ + v^ + w 2) 

Ai 2V > 

and the pressure is obtained from the thermodynamic relation 

p = pe (y - 1) (12) 

where y is the ratio of static enthalpy h to internal energy e (y = h/e) with the heat 
of formation of the gas adjusted so that h = e = 0 at a temperature of absolute zero. It 
should be noted that equation (12) is an exact relation that holds for any gas and is dis- 
cussed in reference 13. Also, it is shown in that reference that the speed of sound for an 
equilibrium gas is related to e and y(p,e) by the equation 
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a2 = e 


(13) 


(9 - i) 


y + 


by 


9 log e e. 




8y 


loge PJ 


For a calorically perfect gas, y is equal to the ratio of specific heats, but these quanti- 
ties differ for real gases. A curve fit for the thermodynamic function y for equilibrium 
air in terms of the density p and internal energy e is given in reference 13 and is 
also presented in appendix A. Expressions for the partial derivatives of y in equa- 
tion (13) are obtained by differentiating the expressions for y given in this appendix. 

The Von Neumann conditions provide a means of estimating the permissible values 
of At and e which can be used at a given time step. It is shown in appendix B that the 
damping coefficient e must satisfy the inequalities 


0 ^ e 


_ J_ 

max ” 24 


and that the time step At must satisfy the inequality 


At S . 


f u (e) min(A AX, 6 AY,r A4>) 


(14) 


/3(v + a)|l + 1 Q(Q + ^4 + Q2 


at each point, where V is the magnitude of the total velocity, e = e /e max , f u (e) is a 
function which is derived and discussed in appendix B, and Q is given by the equa- 
tion (B23). For each time cycle, the value of At is taken to be the smallest value of 
the right side of inequality (14) which occurs at any point in the flow field at the previous 
time step. 


Procedure at axis of coordinate system.- The governing partial differential equa- 
tions (9) contain indeterminate forms on the axis of symmetry, and the velocity compo- 
nents u and w are multivalued with respect to the angle $ there. In this paper, the 
solution for X = 0 and $ = 0 is determined with finite-difference equations, and the 
solution at the axis for $ not zero is determined algebraically. Indeterminate forms 
appear in the governing equations at the axis where X = 0. When these terms are evalu- 
ated with l'Hospital's rule, equations (9) are written as 


9Ai 

Tt 




a 2 Di — 
+ - + Ei 

8X 8 $ 
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where 


A 


Jpl 

II 

CO 

(i = 1,3,5) 

Ej = Kpv 

B2 = p + 2pu^ 


E2 = 2Kpuv 




( 9 a\ 

Dj = pw 


e 3 = -k 

p + plu^ - V") 
J 

D2 = puw 
D3 = pvw 
Dg = pwH 


E 5 = KpvH 

J 


The cross-flow momentum equation is not needed because it is known from symmetry 
that w = 0 in the plane <i> = 0. 

The flow properties for X = 0 and $ / 0 are determined from the kinematic 
relations 


P( t jO,Y ,<3?) = p(r,0,Y,0) 
u(T,0,Y,«f>) = u(t,0,Y,0)cos $ 
v( T ,0,Y,<&) = v(t,0,Y,0) \ 

w(r,0,Y,<3>) = -w(T,0,Y,0)sin $ 
e(r,0,Y,$) = e(r,0,Y,0) 


(16) 


Procedure at downstream boundary .- The region of computation terminates at the 
downstream boundary X = X max where the flow must be supersonic. The procedure 
which is used at points on this boundary is the same as that used at general points except 
that the derivatives 8B i/^X in equations (9) are evaluated with backward differences of 
the form 

~ Bi (r ,X max , Y , =^|®i( T - X max,Y^) - Bi(r;x max -AX,Y,.3>)] 
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Procedure where curvature of reference surface is discontinuous.- A special pro- 
cedure is used for values of X where the curvature of the reference surface changes. 

As shown in figure 2, these are the locations of the junctions of the segments of constant 
curvature comprising the generator of the reference surface. The procedure consists of 
extending the grid on one side of the line of discontinuity or the other by one mesh 
spacing, as shown in figure 4, and performing the calculation in a mesh block with con- 
stant reference-surface curvature. The choice of which grid to extend is determined 
from the arc lengths along the line Y = Constant and $ = Constant to the neighboring 
mesh points; the grid associated with the shortest arc is the one which is used. Thus, 
the polar grid is used to make calculations for point A since the length of arc AB is 
less than that of arc AC, and the alternate grid is used at point D since the length of 
arc ED exceeds that of arc DF. 

Calculations are made at points A and D with the appropriate curvature K 
and equations (9), (10), and (11) as in the general case. Consider point A. The partial 
derivatives with respect to Y and 4> are formed in the usual manner. The partial 
derivatives with respect to X at A are determined from the differences between the 
function Bj at points C' and B. The functions Bj at point C f are determined by 
interpolating the values at points B and C. The interpolation technique is illustrated 
in figure 5. First, the components of the velocity in the plane <3? = Constant at points B 
and C are expressed in terms of the components Ug,Vg and Uq,v^ which are paral- 
lel to the components UQ t ,VQ f . Then the functions Bj at points B and C are deter- 
mined in terms of u,v and the other flow properties. Finally, the values of the func- 
tions at C' are determined by interpolating between the points B and C. It should 
be noted that interpolating linearly between points A and C leads to numerical 
instability. 


Calculation of Flow at Body Points With Continuous Slope and Curvature 

The procedure which is used at body surface points where the slope and curvature 
are continuous is similar to that used at points in the shock layer in that calculations are 
made with equations (9), (10), and (11). The major differences are that the derivatives 
with respect to Y are evaluated with one-sided second-order finite-difference expres- 
sions of the form 




-Ci(r,X,2 AY,<3>) + 4Ci(r,X,AY,$) - 3Ci(r 


,X,0,4>) 


and that the damping terms 


(AY) 4 


8 4 Ai 


8Y- 


in equations (11) are replaced by terms of the 


form 
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2 

(AY) 2 JL- A i (r-AT,X,0,$) = Ai(T-AT,X,2AY,<f>) - 2Ai(r-AT,X,AY,4>) 


BY* 


+ Aj(t-At,X,0,4>) 


where Y = 0 at the surface. At points adjacent to sharp shoulders and at points adja- 
cent to curvature discontinuities where the derivatives of the quantities Bj undergo 

d A.* 

large changes, the damping terms (AX)^ in equations (11) are replaced with terms 


9X^ 


of the form 


-(AX) 2 -ii- Ai(T-AT,X,0,S) 
8X 2 


— Aj 


(t-At,X+AX, 0,$) - 2Ai(T -A t,X,0,«E>) 
+ Ai(T-AT,X-AX,0,<f>)J 


The nature of flow at curvature discontinuities is discussed subsequently. The proce- 
dures which are used at the axis of symmetry and the downstream boundary are the same 
as those used in the shock layer. 

Special considerations must be given to the calculation of flow at the surface when 
there is a discontinuity in either the surface slope or the surface curvature. These con- 
siderations are discussed subsequently. 


Analysis of Flow at Sharp Shoulder 

Consider the flow past an axisymmetric body with a sharp shoulder at angle of 
attack. At the surface there are two nonzero velocity components, u and w, which are 
tangent to the surface and in the radial and cross-flow directions, respectively. At the 
shoulder the u-component is normal to the edge, and the w-component is tangent to it. 

B ehavior of cross flow at shoulder .- Consider an axisymmetric body with a sharp 
shoulder such as that shown in figure 2 or 4. It can be shown that the derivative 8w/9x 
is bounded at a sharp shoulder for steady flow. The steady cross-flow momentum equa- 
tion is obtained from equation (2d) as 


9w 

dx 


1 

r 



P d(pj 


+ w cos e 


(17) 


The quantities w, 1/p, and 1/r are known to be finite at the shoulder, and there is no 
reason for the derivatives with respect to cp to be unbounded. It is shown subsequently 
that the u-component of velocity must be at least as large as the sonic velocity at the 

19 



shoulder so that the quantity 1/u is finite. Therefore, it can be concluded that the 
derivative of w at the shoulder is finite. In this paper it is assumed that this charac- 
teristic applies to unsteady as well as steady flow. 

Sonic condition for n onzero cross flow.- It is well known that for axisymmetric and 
plane flow the surface velocity must either be supersonic upon arrival at the shoulder or 
become sonic there. It can be shown with equations (2) and (4b) that for three-dimensional 
flow it is the surface component of velocity u, normal to the edge, rather than the total 

surface velocity yu2 + w2 which must satisfy one or the other of these conditions. It 

can be seen from figure 2 that in the vicinity of the shoulder ^x c m j n = x = x c max j the 

curvature of the reference surface K, the scale factor X, and the perpendicular distance 
from the axis r are written as 


K = 




r = r c + y sin 8 


respectively, where y = y - y^. By using these relations and equation (1), the equations 
of motion can be written as 


3P , u 9P , v 8p , w dp 
3t y 88 9y r c + y sin 8 88 


+ P 


^ 3u + ^V + 1 8w 

y 88 8y r c + y sin 8 8cp 


cos 8 

r c + y sin 8 



sin 8 


r c + y sin 8 


= 0 


3u u 3u 3u w 9u 1 9p uv 

— + — + V — + — + — — + — 

9t y 3 8 9y r c + y sin 8 8<p py 9x y 


cos 8 

r c + y sin 8 


w 2 


= 0 


i^ + HiX + v — + w j^v + i . sin 8 w 2 = 0 

9t y 88 8y r c + y sin 8 8 cp p 9y y r c + y sin 8 


9w u 8w v 9w + w 9w 1 9p_ 

9t y 3 8 3y r c + y sin 8 8<p p(? c + y sin dj 3 cp 

u cos 8 + v sin 8 w _ q 
r c + y sin 8 
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3p , u Bp , „ 9p , w 3p | a 2 p 1 9 u , 3v , 1 9w 

9t y 90 9y r c +ysin 0 d<p y 90 9y r c + y sin 9 dcp 


cos 9 1 sin 6 V n 

+ u + — + v =0 

r c + y sin 6 \y r c + y sin 9 


Let the flow properties in the vicinity of the corner be represented by functions of the 
form 

v( t,9,0) 

F{t,9,y,cp) = F c (t,0,^) + y F 0 (t,y ,9,cp) (18) 

where F 0 is bounded for all values of y and the exponent v is greater than or equal 
to zero. When these functions are substituted into the equations of motion, the resulting 
equations are multiplied by y, and the limit is taken as y approaches zero, the terms 
involving v and F 0 are eliminated, and the following equations are obtained: 


9p c /9u r \ 

U c T9~ + Pc W + v c = 0 


/9u \ 1 9p 

i c — + v c +- -= o 

C \9 9 C / p c 99 


- u r = 0 


u c — ^ = 0 
c 9 9 


9 p c 2 r u c > 

i c + a^Pp + Vp 

c 9 9 c c \9 9 i 


Equations (19a), (19b), (19c), and (19e) are the Prandtl-Meyer equations where, by 
assumption, p c , p c , u c , v c , and a c are functions of t and cp as well as 9. 
From these equations, it can be shown that the following condition must hold: 

- 4 ) - 0 


I 


Thus, the u-component of velocity either must be sonic for x c min = x = x c max or it 

must be supersonic so that the constant-density or wedge-flow solution applies. The 
solution to equations (19a), (19b), (19c), and (19e) for a perfect gas can be found in 
appendix B of reference 9. In order to use this solution for angle-of-attack flow fields, 
the quantity in reference 9 should be interpreted as follows: 


x=x c,min 

Equation (19d) simply states that the cross-flow velocity component is not a function of 6 
at the shoulder. 

Differentials of pu, p + pu 2 , a nd puH.- Consider isentropic, isoenergetic flow on 
a surface approaching a sharp shoulder. The differential energy equation can be written 
as 

dH = u du + w dw + dh = 0 


where 


1 a 2 

dh = - dp = - 2 - dp 

P P 

Thus, the differential du can be written as 

du = - - dw - — dp 
u pu 

It follows that the differentials of the quantities Bj, B 2 , and B 5 can be written as 
d(pu) = i(u 2 - a 2 ) dp - p ^ dw 

d^p + pu 2 j = ^u 2 - a 2 ^dp - 2 pw dw > ( 20 ) 

d(puH) = -jj(u 2 - a. 2 ) dp - p — dw 
u \ u j 

At a sharp sonic shoulder the quantity (u 2 - a2)dp in these equations is an indeterminate 
form and can be evaluated by expanding the thermodynamic properties and the velocity 
components in terms of the distance from the corner ? = x c rn j[ n - x. It follows from 
equation (17) that the cross-flow velocity component w varies linearly with £. The 
expressions for the density p, the speed of sound a, and the velocity component u are 
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( 21 ) 


p = p*+p iC w 

a = a* + aj? w ) 

u = a* - uj? 40 

J 

where w is a number between 0 and 1 and the coefficients pj, aj, arid uj are posi- 
tive. It has been shown by Guderley (ref. 14) for irrotational flow that co has a value 
of 2/5. Friedman (ref. 15) has shown that this value of co also applies for rotational 
flow. It can be shown with equations (21) that the indeterminate form can be written to 
lowest order in £ as 

^u 2 - a 2 )dp = -2wa*p 1 ^u 1 + 

Thus, the indeterminate form is, in fact, infinite at the shoulder since the exponent 
2co - 1 has a value of -1/5. As a result, the differentials of Bj, B2, and B5 in 
equations (20) are also infinite at the shoulder. Terms of order 3co - 1 and higher do 
not contribute to the values of the differentials at the shoulder since 3a> - 1 has a value 
of +1/5. It can be seen from equations (21) that the differentials of the thermodynamic 
properties and the u-component of velocity vary as = £-3/5 and, hence, have 

stronger singularities than the differentials of Bj, B2, and B5. 


Methods of Calculation at Sharp Shoulders 

Two computational procedures are used at sharp shoulders. One procedure is 
employed when the flow at the shoulder is supersonic and hence nonsingular, and the 
other is used when the flow is sonic and hence singular. The choice of procedure is 
made on the basis of the value of the ratio u/a at the previous time step. 


Procedure used at sharp supersonic shoulder .- If the u-component of velocity at 
the shoulder at the previous time step is supersonic, equations (9), (10), and (11) are 
used as in the general case and the derivatives dB\^BX in equations (9) are evaluated 
with backward-difference expressions. The value of the curvature used in equations (9) 
is that of the segment of the reference line located upstream of the shoulder. Thus, only 
the body geometry and the flow upstream of the shoulder influence the solution at the 

94Ai 


shoulder, as should be the case for supersonic flow. The damping terms (AX) 4 


0X 4 


in 


equations (11) are replaced with second-order terms of the form 
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( 22 ) 


(AX) 2 


Ai(T-AT,X,0,*) 

ax z 


= A i (r-Ar,X-2AX,0,*) 


2Ai(r-Ar,X-AX,0,$) 


+ Ai(r-Ar,X,0,#) 


Procedure used at sharp sonic shoulder .- The procedure described for use at a 
sharp supersonic shoulder was tried also at a sharp sonic shoulder. As is shown in the 
section on "Results and Discussion," the solution obtained at the shoulder in this manner 
is subsonic and does not exhibit singular behavior. It is apparent that special considera- 
tion must be given to the calculation of flow at a sharp sonic shoulder. 

The difficulty with using equations (9) for i = 1, 2, and 5 at a sharp sonic shoulder 
may be due to the fact that each of these equations has two derivatives which become 
unbounded as a sharp sonic shoulder is approached along the upstream surface in such a 
manner that the sum of the derivatives is always finite. It has already been shown that 
the derivatives SBi/sX for i = 1, 2, and 5 become infinite as a sharp sonic shoulder 
is approached along the upstream surface. It can be shown that the same behavior is 
exhibited by the derivatives 9Cj/9Y for i = 1, 2, and 5 since the derivative 9v/8Y 
at the surface becomes infinite as a sonic shoulder is approached from upstream (see 
ref. 15 for example) and since the scale factor X does not vanish on the upstream sur- 
face. Instead of using equations (9) for i = 1, 2, and 5 at sharp sonic shoulders, a set 
of approximate time-dependent equations are employed. 

A set of exact equations in terms of derivatives with respect to X and $ which 
is applicable at the body surface for steady flow can be obtained with equations (20), (17), 
and (8). For example, the first of this set of equations is written as 



This equation is one form of the steady continuity equation, and the second and third 
equations in the set are forms of the steady x-momentum and the energy equations, 
respectively. Therefore, the new equations must be the same as equations (9) for 
i = 1, 2, and 5 in the time-asymptotic limit. 

Two basic assumptions are made in order to obtain the approximate time -dependent 
equations which are used at sonic shoulders in this paper. First, it is assumed that the 
new set of equations is applicable to transient as well as steady flow when the appropriate 
time derivatives are added. Second, it is assumed that the terms proportional to the 

quantity (u 2 - a?)— can be neglected. The approximate time-dependent continuity, 
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x-momentum, and energy equations which are used at sonic shoulders are written as 


8 (pjQ _ f 8 (P u) _ PW fl L 8 w [ 1 9p\ 
St 8X ru [u\ 8$ p ?,§) 


+ w cos 6 




8 (puX) _ J 

fs(p + pu 2 ) 2pw 

El 


8t 1 

1 8X r 

g 1 

y 8<f> p 8 <f>) 


+ w cos 6 




a [(pH - p)x] 

8r 


fa(puH) pwH 

1, 

( w ®z + I®e} 

l 8X ur 

u 

( 8$ P 8$/ 


+ w cos 9 


( 23 ) 


respectively. The value of the reference-surface curvature K which is used in the 
expression for X is that associated with the body surface upstream of the shoulder. 

The previous version of this method presented in reference 7 differs from the present 
version in that equation (17) is not used to replace the derivative 8w/8X which other- 
wise appears in equations (23). 

The effect of neglecting the terms proportional to the indeterminate form may be 
to change the order of the singularity. If the thermodynamic properties and the 
u-component of velocity near a sonic shoulder vary as £l/2+c*-> (3j is a small non- 
negative number) rather than as £2/5 as i S physically correct, the quantity (u2 - a 2 J— 

would vary as rather than as £"V5 and hence would vanish at the shoulder. 

The numerical procedure which is used at a sharp sonic shoulder consists of inte- 
grating equations (23) and the standard cross-flow momentum equation (i = 4 in eqs. (9)) 
with the two-step numerical technique given in equations (10) and (11). The cross-flow 
momentum equation is used in its standard form since w is not singular at the shoulder. 
No damping terms are used in the integration of equations (23). 


Calculation of Flow at Curvature Discontinuities 

Singular behavior can also occur at body points where the curvature of the surface 
is discontinuous but the slope is continuous. These are the junctions between the seg- 
ments of constant curvature which compose the generator of the surface. It should be 
remembered that the only bodies with discontinuities in slope or curvature which are 
treated with the present method are axisymmetric bodies. Two examples are shown in 
figure 6. The singular behavior consists of discontinuities in the partial derivatives of 
p, p, and u with respect to X and it can occur for both subsonic and supersonic flow. 
In general, the streamwise derivatives of the flow properties are discontinuous where the 
surface curvature is discontinuous. The changes in the derivatives in subsonic regions 
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where the surface curvature increases abruptly in the direction of flow tend to be much 
smaller than the changes in the derivatives in supersonic regions and in subsonic regions 
where the surface curvature decreases abruptly in the direction of flow. It should be 
noted that the flow tends to compress where the curvature decreases abruptly and it tends 
to expand where the curvature increases abruptly. The junctions where the X-derivatives 
tend to undergo large changes are designated by the letter L in figure 6. 

If the flow at the junction is supersonic, it is calculated in the same fashion as at a 
sharp supersonic shoulder. Another procedure is employed when the flow at the junction 
is subsonic. This procedure involves the use of the governing equations in a form which 
does not contain curvature terms explicitly. These equations are 

£P + »(PU) + £ £L + + cos A p u _ o 

ar as 6 ay a$\r ) r 

8(pu) + a(p + pu 2 ) + pu ay ( a /puw\ [ cos e J 2 _ w 2) = 0 

Br Bs 6 BY 8$\ r ) r v j 

8 (pw) + a(puw) + pw av + a /p + pw 2 \ + 2 cos e _ 0 

ar as 6 ay a$\ r / r 

8 (pH - p) + 8(puH) + pH SV_ + a /pwH\ + COS e puiI = 0 

ar as 6 ay a$\ r / r 

where s is the distance along the surface from the axis of the coordinate system in a 
plane of constant $ so that 

8 18 
as ~ A ax 

Since the solution at a subsonic junction is influenced by both the upstream and down- 
stream flows, one-sided difference expressions cannot be used to evaluate the derivatives 
with respect to s. Instead, the upstream and downstream derivatives are evaluated and 
the results are averaged. Let X + and A_ be the values of the scale factor for points 

with values of X greater than and less than that of the junction, respectively. The 

finite-difference expressions for the partial derivatives with respect to s are 
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(i = 1,2,4, 5) 


iL Bi(T,X,0,*) =^-^£-|Bi(T,X + AX,0,4.) - Bi(r,X,0,$)j 


^-pi(T,X,0, 


<f>) - Bi(t,X-AX 




At points where the modified equations are employed, the damping functions which 
are used are of the same form as those used at sharp supersonic shoulders and are given 
by equations (22). 


Calculation of Flow at Bow Shock Wave 


The method of characteristics is used to determine the flow properties at the bow 
shock wave. This method consists of solving simultaneously the Rankine-Hugoniot rela- 
tions for a moving shock and one characteristic compatibility relation. In this paper the 
characteristic compatibility relation is integrated with a finite-difference technique. 

Let the direction cosines of the normal to the shock wave with respect to the x-, y-, 
and ^-directions at the shock point be cos ^ cos /3y, and cos 0^, respectively. 

These quantities can be calculated with the equations 


cos jStjj = 


i 8 y s 

AG 3x 



COS /3y = 


1_ 

G 




cos 


P<p 


J_!ls 
rG 3 cp 



(24) 


where 


G = 





+ 1 


( 25 ) 


The unit vector normal to the shock is written as 


The Rankine-Hugoniot relations are written in terms of the thermodynamic vari- 
ables p, p, and y = p the velocity component normal to the shock V, the component of 

the shock velocity in the y-direction 6, the quantity /3y, the angle a between the nor- 
mal to the shock wave and the free-stream direction, and the magnitude of the free- 
stream velocity V M as 


p(v - 6 sec |Sy )j = pjvoo cos a - 6 sec 0yj 
p + p(V - 6 sec /3yJ 2 = p ro + pJV" cos a - 6 sec /3yj 2 

+ 1(V . 5 sec % ) 2 =^7 T ^ cos a - 6 sec 


(27) 


The unsubscripted quantities on the left side of equations (27) are evaluated immediately 
behind the bow shock wave, whereas the quantities subscripted with °° on the right side 
of these equations are evaluated in the free stream. It should be noted that the normal 
velocity V is related to the velocity components u, v, and w by the equation 

V = e’g • V = u cos /3x + v cos /3y + w cos (28) 

where the total velocity vector V is expressed as 

V = e^u + e^v + e^w (29) 

Let the direction cosines of the free-stream direction with respect to the x-, y-, 
and (^-directions of a point be cos a x , cos oiy, and cos a^, respectively. Expressions 
for these quantities are 


cos a x = sin 8 cos a + cos 6 cos cp sin a 
cos Ofy = -cos 8 cos a + sin 8 cos cp sin a > 
cos = -sin cp sin a ^ 


(30) 


where a is the angle of attack. The unit vector in the free-stream direction is 

e— = e x cos a x + cTy cos ay + e ^ cos a^ (31) 

Voo " ' ^ 
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The cosine of the angle 


cr is 


cos a = 



• es = cos a x cos fa + cos ay cos /3y + cos a ^ cos 


(32) 


The characteristic compatibility relation which is used in this paper is derived in 
appendix C. It is written in the form 



' pa2 ^ IF + r 13 2 ) - fK av ' u2 C ° S % + UT cos 


pa 

r 


^av - w^ cos jSy + vw cos jS^sin 6 + (au - w^ cos fa + uw cos p^cos 6 


(33) 

The quantities U x and U<^ in equation (33) are the components of the shock tangential 
velocity component U in the x- and ^-directions, respectively. It can be determined 
from equation (29) and the vector relation 

U = e s x(vx e s ) 


that these quantities and the component Uy satisfy the equations 
U x = u sin? fa - cos fa(v cos /3y + w cos /3^ 

Uy = v sin2/3y - cos fa(n cos fa + w cos /3<^ \ 


(34) 


= w sin 


in^ Pep - COS pylu COS fa + V COS /3y j 


It should be noted that the angles fa, /3y, and (3^ are evaluated at the shock point 
where the flow is being calculated and do not vary with X, Y, or $ when the deriva- 
tives of V, U x , and are being calculated. In other words, the only quantities in 
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equations (28) and (34) which are allowed to vary during the calculation of the derivatives 
are the velocity components u, v, and w. The form of equation (33) which is to be 
used at the axis (X = 0) for # = 0 is 



The solutions at the axis for other values of 4> are obtained from equations (16). 

The present approach differs from the standard time -dependent method-of- 
characteristics approach at shock waves in that the characteristic compatibility relation 
given by equation (33) or (35) is integrated with finite-difference expressions. Central- 
difference formulas are used to evaluate the partial derivatives with respect to X and 
<f>, and backward-difference expressions are used to evaluate the derivatives with respect 
to Y. All these difference expressions are correct to second order. It should be noted 
that this modification of the method of characteristics is similar to that used by Masson, 
Taylor, and Foster (ref. 6) at the shock and body. 

A two-step Brailovskaya scheme of the sort used in the shock layer and on the body 
surface is used to perform the integration. For the first step, the direction cosines of 
the normal to the bow shock wave and the free-stream direction and the shock angle a 
are determined by using equations (24), (30), and (32). Then, the right side of equa- 
tion (33) or (35) and the coefficient pa on the left side are evaluated at time t - At. 

Let the right side of equation (33) or (35) be designated as R. The integrated compati- 
bility equation for the first step is written as 

Pi (t,X,1,4>) + p(T- A t ,X, l,$)a(r- A t,X, 1,$) Vj(t ,X, 1,<£) 

= p(t-At,X,1,<&) + p(t-At,X, l,4>)a(T-AT,X, 1,4>)V (t- At,X, 1 ,$) + AtR(t-At,X,1,$) (36) 

The subscript I used in equation (36) denotes the initial solution. This equation is 
solved simultaneously with the RankLne-Hugoniot relations (27) in order to determine 
an initial solution for the thermodynamic properties p, p, and y, the velocity compo- 
nent V, and the shock velocity component 6 at time t. If an equilibrium gas model 
is being treated, it is necessary to perform an iteration in order to obtain a converged 
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solution for the thermodynamic properties. The velocity components u, v, and w are 
determined with the equations 



For the second step, the right side of equation (33) or (35) and the quantity pa are 
evaluated at time r with the initial solution just obtained. The integrated compatibility 
relation is written as 

p(T,X,l,$) + p^T^l^aj^X.l^M^X,!,*) 

= p(T-AT,X,l,$) + P I (T,X,l,$)a I (T,X,l,$)V(T-AT,X,l,$) + ArR I (r,X,l,i>) 

This equation is then solved simultaneously with the Rankine-Hugoniot relations in the 
same manner as for the first step. 

No damping terms are used at the shock wave since they proved to be unnecessary. 


Starting Solution 

It has been shown by a number of investigators that converged solutions for the 
steady supersonic blunt-body problem which have been calculated with time-dependent 
finite-difference methods are independent of the starting solution for all means and pur- 
poses. Thus, the starting solution can be quite approximate. 

In this paper the starting solution is constructed from a specified axisymmetric 
shock wave with its axis in the free-stream direction and a specified surface pressure 
distribution and surface streamline pattern. It is assumed that the shock velocity for the 
initial solution is zero. Therefore, the complete solution at the shock can be determined 
with the Rankine-Hugoniot shock relations, the specified free-stream conditions, and the 
shock geometry. 

The basic surface pressure distribution is the Newtonian one. At sharp corners for 
which the upstream value of the velocity component u is subsonic the pressure is 
adjusted so that u is equal to the local speed of sound. If the Newtonian pressure dis- 
tribution does not indicate a stagnation point on the windward surface of the body (this 
happens, for example, for a flat-face cylinder at angle of attack), the stagnation pressure 
is imposed at the points with the largest Newtonian pressure, and the pressure at the 
other points is increased proportionally. In regions of wind shadow, the pressure is set 
equal to the pressure at the last point on the windward surface in the same a? -plane. 


31 


I III II II 


There is a provision for imposing a lower bound on the pressure at all points in the flow 
field. 

The density at the surface is determined from the pressure distribution and the 
normal shock entropy. The use of normal shock entropy is considered reasonable since 
the results of Shifrin (ref. 16) suggest strongly that the maximum entropy streamline 
wets the surface when the flow field is not symmetric. The magnitude of the velocity V 
at the surface is determined from the pressure, density, and the total enthalpy. The sur- 
face components of velocity are obtained from the total velocity and the Newtonian surface 
streamlines which are used even in regions of wind shadow. 

The flow properties at points between the shock and body are determined by linear 
interpolation. When an equilibrium gas model is to be used, the initial solution is cal- 
culated as a perfect gas. The value of y which is used is the normal shock value. 

RESULTS AND DISCUSSION 

In this section the results of the present method are compared with those of experi- 
ment and other methods for the flow of equilibrium air past the Apollo command module 
at the trim angle of attack and for perfect gas flow past a spherical cap and a spherically 
blunted cone at angle of attack. Both the cap and the cone are terminated with sharp 
shoulders. Results are presented also for perfect gas flow past a prolate spheroid with 
its major axis normal to the flow. In addition, several solutions for the flat-face cylinder 
at zero angle of attack which were obtained by using different finite-difference techniques 
are compared. 


Smooth Bodies 

Apollo command module .- The results of the present method for the bow-shock and 
sonic-line shapes and surface pressure distribution for the Apollo command module are 
compared with those of the inverse method of Webb, Dresser, Adler, and Waiter (ref. 17) 
in figure 7. The angle of attack is 22°, the flight velocity is 6935 m/s (22 754 ft/sec), 
and the altitude is 45.866 km (150 480 ft) in the atmosphere of the earth. For the pres- 
ent calculations, the free-stream properties at this altitude were obtained from refer- 
ence 18 in terms of the geopotential altitude. It should be noted that since the Apollo 
command module is axisymmetric and has a generator composed of segments of constant 
curvature, it is well suited for treatment with the present method. The results of both 
methods are for air in thermodynamic equilibrium. 

The shock- and sonic-line shapes presented in figure 7(a) are those in the plane of 
symmetry. The results of the present method and those of the inverse method (ref. 17) 
for the shock location and the sonic-line location in the windward semiplane (<p = 180°) 
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coincide to within plotting accuracy. It is seen that the results of the two methods differ 
somewhat as to the location of the sonic line in the leeward semiplane ( cp = 0°). 

The results for the pressure distribution for the plane of symmetry (cp = 0° and 
180°) and for the plane normal to the plane of symmetry (cp = 90°) are compared in fig- 
ure 7(b). It is seen that the difference in the sonic-line location in the leeward semiplane 
( cp = 0°) is associated with only a small difference in the surface pressure distribution. 
This may be due to the fact that the value of the quantity y = h/e is close to 1. Since 
the gas is in thermodynamic equilibrium, the value of y varies somewhat over the flow 
field. However, in the subsonic region, the value of y is within 1 percent of the mean 
value of 1. 145. It is seen in the figure that the inverse method predicts that the expan- 
sion at the corner in the plane cp = 90° starts closer to the axis than the present method 
does. It should be noted that difficulties are often experienced with the inverse method in 
the vicinity of shoulders. 

The calculation presented in figure 7 for the Apollo command module was performed 
with a grid with seven semiplanes of constant cp so that A4> = 30°, three strips between 
the body and shock so that AY = 1/3, and nine mesh spacings AX along the generator 
of the reference surface. Five of these mesh spacings are on the arc subtending the face 
of the configuration, and the remaining four are on the arc subtending a portion of the 
shoulder. The time required to calculate the results of the present method shown in fig- 
ure 7 on the CDC 6600 computer was 10 minutes. However, the results obtained after 
5 minutes of calculation did not differ appreciably from the results presented in the fig- 
ure. The inverse solution of reference 17 required many hours of computer time on the 
IBM 7094. A second calculation was made with the present method with a more refined 
grid. The results of this calculation did not differ appreciably from those presented in 
figure 7. 

It was noted that the numerical solution is extremely sensitive to errors incurred 
by using too few ^-planes. An attempt was made to calculate the flow field about the 
Apollo command module by using the same grid as that just discussed except that five 
rather than seven planes of constant cp were used. This attempt was unsuccessful. 

The results obtained indicated that a pseudoseparation occurred on the shoulder on the 
leeward side of the body. First the cross-flow component of velocity w and then the 
u-component reversed directions near the shoulder. The total enthalpy increased in the 
region of reverse flow. Eventually, the calculation became unbounded. 

Prola te spheroid .- In order to demonstrate the ability to compute flow about non- 
axisymmetric bodies, two calculations with different mesh spacings were made for a 
prolate spheroid with its axis perpendicular to the flow. It should be noted that a prolate 
spheroid is a body of revolution with an ellipse for a generator and with the axis of revo- 
lution passing through the foci. For the first calculation, seven planes of constant cp 



and three strips between the body and shock were used so that A<h = 30° and AY = 1/3. 
Seven mesh spacings AX were used along the generator of the reference surface. For 
the second calculation, the mesh spacings A<£ = 22.5° and AY = 1/5 were used. 

Along the generator of the reference surface, 14 mesh spacings AX were employed. 

The bow-shock-wave and sonic-line shapes in the planes cp = 0° and <p = 90°, which 
are normal to and contain the major axis, respectively, are shown in figure 8(a). The 
surface pressure distributions for these planes are shown in figure 8(b). Since the 
results of the two calculations are in good agreement, it can be concluded that conver- 
gence has been attained. Also shown in figure 8(b) is the Newtonian distribution. It is 
seen that the present results and the Newtonian results are in good agreement for the 
plane cp = 0° but differ somewhat for cp = 90°. 

Bodies With Sharp Shoulders 

Flat-face cylinder . - In the Analysis section of this paper, it was indicated that a 
special procedure was needed in order to use finite-difference equations at sharp sonic 
shoulders. In order to demonstrate the need for the refinement, calculations were made 
for the flow past a flat-face cylinder at zero angle of attack with Moo = 2.81 and y = 1.4 
with both the modified and the unmodified procedures. The results for the surface pres- 
sure distribution and the pressure profile across the shock layer along a line normal to 
the face at the shoulder are presented in figure 9. The circles depict the results of the 
present method, which uses equations (23) at the shoulder when the upstream flow is sub- 
sonic, and the squares depict the results of a variation of the present method which uses 
equations (9) at the shoulder. The results depicted by the triangles were taken from ref- 
erence 9 where the method of characteristics was used at the body surface. The experi- 
mental results of Kendall (ref. 19) are depicted by a solid line. 

Since the present method employs equations (9) at all points except the shoulder, it 
might seem appropriate to use these equations at the shoulder also. However, it is seen 
from the figure that such a procedure does not predict the solution correctly at the shoul- 
der. Both the present method and that of reference 9 are in fair agreement with experi- 
ment and each other, and the singular character of the flow at the shoulder is predicted 
correctly. It should be noted that the error in the results obtained with equations (9) at 
the shoulder is restricted to the neighborhood of the shoulder. In fact, it was found that 
the results of both the present method and its variation for the shock-wave and sonic-line 
locations agree with those of experiment and method of reference 9. 

Spherical cap .- In figure 10, the results of the present method for a spherical cap 
with a sharp shoulder are compared with the experimental results of Stallings and Howell 
(ref. 20). The angle of attack is 15°, the free-stream Mach number is 2.49, and y is 
1.4. The radius of curvature of the cap is ^r^. The shock-wave and sonic-line 
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locations in the plane of symmetry are shown in figure 10(a), the pressure distributions 
for the plane of symmetry (cp = 0° and cp = 180°) and the plane normal to the plane of 
symmetry (cp = 90°) are given in figure 10(b), and the distributions of the Mach number 
and the ratio u/a along the shoulder edge are given in figure 10(c). The mesh spacings 
A<& and AY have values of 22.5° and 1/3, respectively. There were 15 mesh spacings 
AX along the generator of the reference surface. 

The interest in calculating this particular case stems from the fact that experimen- 
tal results show the flow at the shoulder to be sonic on the windward side of the body and 
supersonic on the leeward side. Therefore, this case provides a good test of the ability 
of the present method to calculate both sonic and supersonic shoulder flow. It is pointed 
out in the Analysis section that different computational procedures are used at sharp 
shoulders, depending upon whether the upstream flow is subsonic or supersonic. When 
the Mach number of the approaching flow is less than 1, equations (23) are used; when it 
is greater than 1, equations (9) are employed. It can be seen from the results presented 
in figure 10 that the present method predicts the character of the flow in the vicinity of 
the shoulder correctly. 

S pherically blunted cone with sharp shoulder .- The results for the flow about a 
spherically blunted cone with a sharp shoulder are compared with the experimental 
results of Stallings and Tudor (ref. 21) in figure 11. The semiapex angle of the cone is 
60°, the ratio of nose radius to base radius is 0.25, the Mach number is 4.63, and the test 
gas is air. The values of angle of attack which are treated are 10° and 20°. The grid 
which is used for the a = 10° calculations has 25 mesh spacings AX along the gen- 
erator of the reference surface, four mesh spacings AY between the body and shock, 
and six mesh spacings A4> between cp = 0° and cp = 180°. The grid used for a = 20° 
has the same number of mesh spacings AX and AY and eight mesh spacings A4>. 
Converged results were obtained for each of these cases in 300 time steps which corre- 
sponds to about 30 minutes of computer time. 

The pressure distributions in the plane of symmetry are shown in figure 11(a). In 
addition to the results of the present method and experiment, the approximate solution 
obtained with a time-dependent, finite-difference method by Cohen, Foster, and Dowty 
(ref. 5) for a 10° angle of attack is also presented. This method, which was discussed 
briefly in the Introduction, is basically the refined Godunov method of Masson, Taylor, 
and Foster (ref. 6). The flow is calculated in the plane of symmetry only, and the cross 
flow is assumed to vary as sin cp. Two solutions for different values of the cross-flow 
weighting function are presented in reference 5. The solution presented in this paper is 
that designated as the ".433 crossflow approximation" in reference 5. 

Several observations can be made concerning the results presented in figure 11(a). 
For the 10° angle of attack, both the present method and the approximate method of 
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reference 5 show good agreement with experiment. Both methods indicate that the flow 
stagnates at a point on the spherical cap in agreement with experiment. It is seen that the 
present results indicate that the pressure on the cap is about 2 percent less than the exact 
stagnation value p t . It should be noted that for the present method, four mesh spacings 
AX were used between the geometrical stagnation point and the edge of the spherical cap. 
It is believed that better agreement with experiment would have been achieved had more 
mesh spacings been used on the cap. Six mesh spacings were used on the cap for the cal- 
culation of reference 5, and excellent agreement with experiment was achieved. It is seen 
that a discontinuity in the pressure gradient at the junction of the cap and cone in the plane 
of symmetry is predicted by both methods for both the windward and leeward sides. 

The results of the present method for the 20° angle of attack are in good agreement 
with experiment except near the shoulder on the leeward side where an irregularity in the 
present results for the pressure distribution occurs. This irregularity occurs concur- 
rently with a nonphysical variation of the surface entropy of 4 percent near the corner. 

The irregularity does not affect the Mach number distribution adversely as it does the 
pressure distribution. The computed Mach number at the leeward corner is 0.969, which 
is reasonably close to 1, and the Mach number computed at grid points in the leeward 
plane upstream of the shoulder decreases monotonically as it should but remains close to 
1 for some distance. The present results show that the stagnation point is located on the 
conical surface in the windward plane in agreement with experiment. 

It was stated in reference 5 that a calculation was made for a 20° angle of attack but 
that the results did not agree with experiment. It was reported that the stagnation point 
was located on the spherical cap in the windward plane rather than on the cone. The 
results of the present method give an indication as to why the approximate method of ref- 
erence 5 was successful for a 10° angle of attack and unsuccessful for 20°. The periph- 
eral distribution of the cross flow along the edge and along the junction of the cap and 
cone are shown for 10° and 20° angles of attack in figures 11(b) and 11(c), respectively. 
Also shown in these figures are the distributions of the form 

W = W 90° Sin 9 

which is the form assumed in reference 5. It is seen that for a 10° angle of attack, the 
computed cross flow varies as sin <p both at the edge and at the junction. For a 20° 
angle of attack, the cross flow at the edge varies as sin <p, but at the junction in the nose 
region the computed cross-flow solution and the trigonometric function differ slightly. 
Apparently the solutions for the various flow properties in the stagnation region are 
strongly coupled so that even a small deviation in the cross-flow velocity distribution has 
a marked effect on the other properties. 
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The peripheral pressure distributions for three values of s/r^ are shown for 10° 
and 20° angles of attack in figures 11(d) and 11(e), respectively. The values of s/t ^ 
are 0.10, 0.45, and 0.90 and were chosen because there are both computational points and 
experimental orifices very close to these locations. The first ring (s/r^ = 0.10) is 
located on the spherical cap near the junction with the cone, the second ring (s/r^ = 0.45) 
is located on the cone and is closer to the junction than the edge, and the third ring 
(s/r^ = 0.90) is located on the cone near the edge. It is seen in figures 11(d) and 11(e) 
that the results for both angles of attack are in fair agreement with experiment but that 
the present and experimental results for the 10° angle of attack are in closer agreement 
than those for 20°. It should be pointed out that at the location s/r^ = 0.10 where the 
comparison is made on the cap, there is more difference between theory and experiment 
than at the other computational points on the cap, as can be seen in figure 11(a). 

In figure 11(d), results are presented for two calculations for the 10° angle of 
attack. For one calculation the mesh spacing A# has a value of 30°, whereas for the 
other calculation the value of A<f> is 45°. The number and size of the mesh spacings 
AX and AY are the same for both calculations. The results of the calculation with 
A4> = 30° are the ones presented in the rest of figure 11. A comparison of the results 
of the two calculations presented in figure 11(d) shows that the solution is converged near 
the edge (s/r^ = 0.90) for all values of cp. It is also converged on the leeward side for 
all values of s. However, it is seen that a refining of the grid does improve the results 
at other points on the surface. 

It was shown in the Analysis section that the cross flow should not exhibit singular 
behavior at sonic shoulders although the other flow properties do show singular behavior 
at these shoulders. In order to show that the numerical solution does, in fact, behave in 
this manner, the results at 10° and 20° angles of attack for the pressure, velocity com- 
ponent u, and cross-flow component w in the plane cp = 90° are shown in fig- 
ures 11(f), 11(g), and 11(h), respectively. It is seen that both the pressure and the veloc- 
ity component u have steep gradients at the shoulder whereas the cross flow does not. 

CONCLUDING REMARKS 

A two-step, time-dependent method of second-order accuracy for calculating super- 
sonic flow about nonaxisymmetric blunt bodies and axisymmetric blunt bodies with sharp 
shoulders is presented. The bow shock wave is treated as a discontinuity, and a modifi- 
cation to the time-dependent method of characteristics is used to determine the solution 
at the shock. Finite-difference techniques are used to calculate the solution at the body 
surface and in the shock layer between the shock and surface. An approximate finite- 
difference method for use at sharp sonic shoulders is presented. A stability analysis of 
the finite-difference method is given. 


37 


Comparisons of the results of the present method with experiment and the results 
of other methods are made for the flow of equilibrium air past the Apollo command 
module at the trim angle of attack and for perfect gas flow past a spherical cap and a 
spherically blunted cone at angle of attack. Both the cap and the blunted cone are termi- 
nated with sharp shoulders. In general, the agreement with other results is good. 

Results are presented also for perfect gas flow past a prolate spheroid with its major 
axis normal to the flow. 

It is found that time-dependent methods which employ only finite-difference inte- 
gration techniques can be used to calculate numerical solutions at sharp shoulders where 
the flow is singular. In particular, finite -difference techniques can be used to integrate 
the compatibility relation which is employed in conjunction with the Rankine-Hugoniot 
relations in the time -dependent method-of- characteristics solution at the bow shock wave. 
Thus, it is not necessary to construct the characteristics network that would otherwise be 
required. 

The present solutions account for the cross flow at the shoulder. It is shown that 
when there is cross flow at the shoulder, the sonic or supersonic nature of the solution 
depends upon whether the component of velocity tangent to the surface and normal to the 
shoulder edge is subsonic or supersonic upstream of the shoulder. The properties which 
are singular at a sonic shoulder are the component of velocity just discussed and the 
thermodynamic properties; the cross flow is not singular. 

It is demonstrated that numerical solutions for supersonic flow about blunt bodies 
at angle of attack can be obtained in a reasonable amount of computer time. For example, 
the flow field with real-gas effects about the Apollo command module at the trim angle of 
attack was calculated in 10 minutes on the CDC 6600 computer. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., May 25, 1971. 
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APPENDIX A 


CURVE FIT FOR THE THERMODYNAMIC FUNCTION f(p,e) 

FOR EQUILIBRIUM AIR 

The data of Allison (ref. 22) and Browne (ref. 23) have been used to construct this 
curve fit for the function y = h/e for equilibrium air. The independent variables are 
the density p and the internal energy e. The ranges of these variables for which the 
curve fit applies are 

10~ 4 S % 10 —2— = 1500 

Pq RTq 

The heat of formation of the gas must be adjusted so that h = e = 0 at a temperature of 
absolute zero. 

The curve fit is given by the following equations: 


Y = 1.405 




y = 1.5055 - 0.1255 log 1Q - 


' = f 1.6370 - 0.0404 log 10 - fo.2175 - 0.0332 log 1Q £-W 10 
\ p o ) \ p o } RT o 

^0. 1366 - 0.0366 log 1Q - ^0.0833 - 0.0248 log 1Q ^log 1Q ~ 


1 + exp 


-18.3o(log 10 - 0.0320 log 10 f- - 1.943) 


y = [l.5004 - 0.0038 log 10 ^0.1342 - 0.0084 log 10 ^)log 10 

[o.3274 + 0.0091 log 10 ^) - [o.l342 + 0.0016 log 10 ^-)log 10 gfr 


1 + exp 


' 15 - 85 ( log io ^ - °- 0320 log i0 ^ - 2 - 708 ) 


3.255 - 2.278 log 10 ^ 


0.801 < login — = 2.300; log 10 ^ 

10 RT o 10 p o 1 - 0.822 log 10 


3.255 - 2.278 log 10 ^ 


0-801 < log 10 S 2.300; log 10 £ < — 

RT 0 p o 1 * 0.822 log 1 n 

10 RT, 


[2.300 <log 10S |-S 3.176) 


The quantity p Q is the reference density and has a value of 1.292 kg/m3 (2.507 x 10“^ 
slug/ft3). The quantities R and T 0 are the gas constant for air and the freezing 
temperature at 1 atmosphere of pressure; the product RT 0 has a value of 78.40 kJ/kg 
(8.439 x 10 5 ft 2 /sec2). 
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APPENDIX B 


STABILITY ANALYSIS FOR FINITE -DIFFERENCE EQUATIONS 


In this appendix the Von Neumann condition for the linearized form of the present 
finite-difference equations is determined. This condition specifies the upper bound for 
the mesh spacing At above which it may be expected that the magnitude of a given 
infinitesimal error will increase at successive time steps. In general, difference 
schemes which amplify small errors produce solutions which either diverge or contain 
large errors. 

The basic procedure which is used here is similar to that of reference 9 and is 
based on the earlier work of Richtmyer (ref. 24). A technique employed by Van Leer 
(ref. 25) is used to account for the difference in the mesh spacing sizes for the X-, Y-, 
and ^-coordinates. This treatment accounts for the nonorthogonal nature of the X,Y,4> 
coordinate system. A nonphysical dissipation function of fourth order similar to that 
used by Richtmyer and Morton (ref. 12) is employed to avoid neutral stability. 

The partial differential equations (2) and (4b) can be written in matrix form in 
terms of the coordinates r, X, Y, and $ as 


3W t A 8W i B 3W | c 8W | p - Q 


8t X 8X 6 8Y 


0$ 


(Bl) 


where the vectors W and D and the matrices A, B, and C are as follows: 


A 


W = { 


r ~\ 
P 

u 

v 

w 

p 

V J 


pu ££±1 + PV \E + ^^ 


uv — - w 2 cos - 
A. r 


D = ( - u 2 K . w 2 sinj 
' X r 


-\ 


uwi22j + TW £i!!_« 


pua 2 c °s - + pva 2 ^ + si ^ - 


J 


(B2) 


(Equations continued on next page) 


40 



I 


APPENDIX B 



The finite-difference representations of equations (Bl) for the initial and final steps 
of the Brailovskaya scheme are written as 
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W I (r,X,Y,#) = W(t-At,X,Y,$) - < f 


W(t-At,X+AX,Y,3>) - W(t-At,X-AX,Y,$) 


2 AX 


|w(t-At,X,Y+AY, j») - W(T-AT,X,Y-AY,$) [ 



+ D > At 


and 


W(t,X,Y,$) = W(t-At,X,Y,$) - 


A 

A 


W!(t,X+AX,Y,$) - Wj^X-AX.Y,#) 


2 AX 


Wj^X/r+AY,^) - Wj(t,X,Y-AY,<IO 



+ D > At 


- ejw(T-AT,X+2AX,Y,$) - 4W(t-At,X+AX,Y,$) + 6W(t-At,X,Y,4>) 

- 4W(t-At,X-AX,Y,3>) + W(t-At,X-2AX,Y,$) 

+ W(t-At,X,Y+2AY,4>) - 4W(t-At,X,Y+AY,4>) + 6W(t-At,X,Y,4>) 

- 4W(t-At,X,Y-AY,4>) + W(t-At,X,Y-2AY,4>) 

+ W(t-At,X,Y,4?+2A4>) - 4W(t-At,X,Y,4>+A4>) + 6W(t-At,X,Y,$) 

- 4W(t-At,X,Y, 4>-A<E>) + W(T-AT,X,Y,$-2A<i>j] 


(B4) 


(B5) 


The subscript I in equations (B4) and (B5) indicates the initial solution. The matrices 
A, B, and C and the vector D are evaluated at the point (t - At,X,Y,#). 

It is assumed that the solution is of the form 


W(t+At,X±AX,Y±AY,#±A 4>) = W(T,X,Y,$)exp 


5 At ± i(?7 + | + 


'P) 

i 


(B6) 
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where 

n = — £X £ = -^L AY ^=|2La$ (B7) 

L x Ly 

and i is '{-1. The quantities L x , Ly, and are the wavelengths of the error 

solutions in the X-, Y-, and $ -directions, respectively. Define the quantities r j, r2, 
and rg as 


r l = 


At 
X AX 


r 2 = 


At 
6 AY 


At 

r« = 

0 r A<t> 


(B8) 


and assume that At is small so that the term D At can be neglected. By using equa- 
tions (B6) to (B8), equations (B4) and (B5) can be written as 

Wj(t,X,Y,$) = eW (t-At,X,Y,$) = EW(t-At,X,Y,$) (B9) 


and 


W(t,X,Y,$) = gW(T-AT,X,Y,$) = GW(t-At,X,Y,4>) 


(BIO) 


respectively, where the matrices E and G are given by the respective equations 


E = r-jA sin t] + rgB sin £ + rgC sin 


(Bll) 


and 


G = (1 - ef)I - E 2 - iE 


(B12) 


where 


= 4^1 - 


cos r\y + (1 - cos £,) z + (1 - cos ip)‘ 


(B13) 


The quantity I is the identity matrix. The numbers e and g are both defined as 
exp(2 At) although they are not equal since the quantity 2 is not the same for the first- 
and the second-step solutions. The equations (B9) and (BIO) can be written in the forms 


(E - I e)W = 0 


(G - Ig)W = 0 
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For nontrivial solutions of W, the determinant of the coefficients of the components of 
W in these equations must vanish; that is 


E - Ie 


= 0 


G - Ig = 0 


(B14) 


Equations (B14) can be used to determine the eigenvalues e and g of the matrices E 
and G, respectively. 

First consider the matrix E. Let the quantities m, cos pj } cos p 2 , cos pg, 
and U be defined as 

'N 


m = 


-)2 


rj sin 77 - 


r 2 (Y6 T + yjQsin I 


9 o 

+ rgSin^ + 


rg sin \fs 


r 2 (^ + y b )sin £ 


r j sin rf - 




sin | 


cos M 1 = 


m 


cos 


M2 = 


r 2 sin £ 
m 


COS (JL 3 = ■ 


sin i// 


r 2( Y ^ + y b )sin £ 


m 


U = u cos p^ + ^v - Yfijcos p 2 + w cos p 3 


J 


(B15) 


The matrix E can be written in terms of these quantities as 


U p COS jUj p cos ju. 2 

0 U 0 


E = 


0 


0 


U 


p cos pg 0 

cos u< 

0 

P 

COS tin 

0 

p 


00 0 u 

0 pa^cos pj pa^cos p 2 pa2cos pg 


COS Pg 


p 

U 


(B16) 
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By using the first of equations (B14) and equation (B16), the eigenvalues of the matrix E 
are 


e = me = mU,mU,mU,m(U + a),m(U - a) 


(B17) 


It should be noted that 


-max 


It can be shown from equations (B12) and (B14) that the eigenvalue of the matrix G 
can be expressed as 


g = 1 - ef - e 2 


le 


The stability condition which is used in this paper is that the magnitude of g 
should not exceed 1. This condition is expressed as 


g ^ = (l - e 2 - ef)^ + e 2 s i 


(B18) 


Consider what happens if e = 0, and note from equation (B13) that 0 = f = 48. It is seen 
from these conditions and the inequality (B18) that e must satisfy the inequalities 


0 = e = e max - 24 

The inequality (B18) can be rewritten as 


(B19) 


MM-f 


+ ef = e2s(i-ef] + 


- + ef = F 




u 


(B20) 


where and Fu are the lower and upper bounds of e 2 . It can be shown that for the 
values of e in inequality (B19) the lower bound has a maximum value of zero. Thus the 
inequality can be rewritten as 


0 = m^e 2 = F 


u 


(B21) 
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From the first of equations (B 15 ) it is seen that m must satisfy the inequality 


m 2 g lr jjsin 77 J + 


sin | 


Y6 + vL 


y b . _2 . 


+ + Irg sin ip + ■ 


r 2 

sin £ 

Y6 + y b ) 

r / 


= m^r 2 


max 


sin 2 77 + sin 2 | + sin?\f/ 


where 


m 2 = ^cos + Pj cos v^j + cos v \ + ( cos v 3 + p 3 cos v 2 j 


COS J'j = 


r j sin 77 


^r jSin 2 77 + r 2 sin 2 £ + r 2 sin 2 i// 
sin d 


COS 


v 2 = 


^rjSin 2 77 + r^sin 2 ! + r^sin 2 !// 
sin i^l 


COS fo = 


\rjSin2?7 + r^sin 2 ^ + r^sin 2 !// 


P l = 


P 3 ^ 


Y6 ’ + y; 


Y6 + y b 


> (B 22 ) 




The quantity r max is simply the largest of r j, r£, and r^. It can be shown that the 
quantity m satisfies the inequality 
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<2. 1.1^2 -2 


1 + COS 2^2) + 2(Pj COS + Pg COS ygjcos ^2 


= K 2 + P 1 - P s) + K P 1 + P 3/ cos 2 ^2 + p! + P 3 sin 2 ^2 

S l( 2 * P 1 + P s) + |H + P 3 I 4 * P 1 * P 3 

= ^ + 1 Q 2 ) + ^ Q\ 4 + Q 2 


where the quantity Q is defined by the equation 


«=H + P 3 = lti^ , + yi ) 2+ ^( YS + ?b 


It can be seen from equation (B13) that the quantity f can be written as 


f = 16[sin^ ^ + sin^ -i + sin^ = 48y2 

V Cl Cl Cl J 


where 0 = x = L Let e be defined as 


^ ~ € / e max ” 


The upper bound F u of e z can be written as 


Fu(x,e) = (I - 26X 2 ) + \g + 2e x 2 


At this point, it should be noted that 


sin2 ^ + sin2 £ + sin2 ^ } = sin^ ^ + sin^ J- + sin^ ^ 

2 2 2 / 2 2 2 


+ 2/sin2 — sin2 — + sin^ — sin^ — + sin^ — sin^ — 

l 2 2 2 2 2 2 


= 3 (sin^ ^ + sin^ ^ + sin^ ^ 

V Cl Ci Cl . 
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It follows that 


i 2 ?7 + sin 2 | + sin 2 i// = 4 ^sin 2 2Z + sin 2 1 + sin 2 - ^sin^ ^ + sin4 1 + sin^ % 12(x - x 2 ) 


Therefore, the inequality (B21) can be rewritten as 



For each value of e between 0 and 1, the function F u (x,e) reaches a minimum value 
for some particular value of x between 1/2 and 1. The minimizing values of e and 
X are related by the equation 

16 X 2 



The least upper bound of F u (x,e) for each value of e is designated as f u (e). This 
function, which varies monotonically from a value of 1 at e = 0 to a value of ][ 2/3 at 
e = 1, is plotted in figure 12. 

In the situation encountered in this paper, the shock velocities are relatively small. 
Thus, the effective upper bound of U is the magnitude of the total velocity V; that is 


U ~ u cos /! i + v cos 112 + w cos [ig = V = \^u 2 + v 2 + w 2 

and the quantity V + a is the effective upper bound of e max = U + a. 

It follows from the analysis presented in this appendix that the Von Neumann con- 
dition for stability for the present finite-difference scheme is expressed by the inequality 

f u (e)min(A AX, 6 AY,r A4>) 

At ? 

l/3(V + a) 
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THE CHARACTERISTIC COMPATIBILITY RELATION 

There are a number of treatments of the multidimensional method of characteris- 
tics. The analysis presented here follows the treatment of Von Mises (ref. 26). 

Consider a multidimensional space. Compatibility relations for this space are 
linear combinations of the equations of motion which can be expressed in terms of total 
derivatives along given lines called a bicharacteristics and partial derivatives in terms 
of some of the independent variables. It should be noted that compatibility relations do 
not exist for all flow fields. However, if one of the independent variables is time and the 
flow field is inviscid, compatibility relations always exist. The total derivatives in these 
relations are with respect to time, and there are no partial derivatives with respect to 
time. 

The space which is used in this paper is four dimensional, and the independent 
variables are time t and the coordinates x, y, and <p. The equations of motion which 
are combined linearly to form compatibility relations along the acoustic bicharacteristics 
are equations (2b), (2c), (2d), and (4b). 

Let the direction cosines of the bicharacteristic in t,x,y,<p space along which a 
compatibility relation applies be ty, h x , hy, and h^. By using the techniques dis- 
cussed in reference 26, it can be shown that there is a family of bicharacteristics through 
each point in t,x,y,<p space, the members of which have direction cosines which satisfy 
the equation 


ht + uhj 


+ v hy + wh^ = ±au 


,2 .2 

h x + hy 


+ h 


<P 


(Cl) 


Let the constants of proportionality for equations (2b), (2c), (2d), and (4b) in the 
compatibility relations be aj, a 2 , a 3 , and a 4 , respectively. It can be shown that aj, 
a 2 , and a 3 are related to a^ by the equations 


a l = - 


pa 2 h x 
h+ + uh Y + vh., + wh 


<P 


a 4 


a 2 = “ 


a 3 = " 


pa^hy 


ty. + uh x + vhy + wh<p 

pa 2 h (p 

h^ + uh x + vhy + wh (p 


a 4 / 


a 4 


(C2) 
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By using equations (Cl) and (C2), the compatibility relations can be written as 


hx /0u u 9u 0u w 0u 1 9p K cos 6 9| 

Tpa — + — — + v — + — — + -- — z- + — uv - - ■ - w* 4 

2 + h 2 + jj2 \ 0t 9x 0y r Bcp pX 0x X r / 


\I4 


t pa- 


ii2 -2 -2 

h x + V\> 


^ + u8v +v 8v + w8v + 18p_K u 2 . sin 6 w 2 
\0t X 0x 0y r 3^ p 9y X r / 


t pa 


V /9w . u 0w , „ 0w , w 0w . 1 0p . cos 6 . sin 6 

— -( 4- v 4 } — 4 UW i 

1 0t X ax dy r 8(p pr d(p r r 


vw 


\ , 2 ,2 ,2 
W*V 


+ )iP + “lE + vl£ + 2iE + p a 2 

|0t X 0x Sy r B(p 


I ^ + + 0w + co^ u + f K + sinj y 

X 0x 0y dcp r \X r / 


= 0 


(C3) 


The quantities 





have all the properties of direction cosines. In this paper these quantities are chosen to 
be the same as the direction cosines of the shock normal, cos cos /3y, and cos (3^, 
so that the total derivatives dp/dt and dV/dt will appear in the compatibility relation 
which is used and dU/dt will not. It has been assumed that the bicharacteristic can be 
approximated by a straight line so that the quantities ty, h x , hy, and h^ can be con- 
sidered constant. With these assumptions, the compatibility relation can be written as 


^ + -i^u + a cos + (v + a cos PyJ-^ + + a cos 



u 0V 0V 
+ + v — 

X 0x 0y 


— — ^ + pa 2 f- — + — + - — \ = --^(av - u 2 cos /3y + uv cos /3J)v - ^ (av - w 2 cos /3y 

r 0 C p \X 9x 9y r Bcp \ \ ** r [V ^ 


+ vw cos /3<^sin 6 + ^au - w 2 cos ^ + uw cos jS^cos 6 


(C4) 
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With the aid of equations (29) and (34), it can be shown that the following relationship 
holds: 


1 8u 8v 1 0W_ cos ftc 8V /j 8V cos Pep 8V 1 8U x 9u y 1 9U <p 

A. 8x 8y r 8<p X 8x y 8y r 9 cp X 8x 9y r 9<p 


Equations (C4) and (C5) can be combined to yield 


M + pa(^ = -pa 2 (l 15* + + I ££<£) - £*L V . u 2 C os /3 V + uv cos 

ydt / + F \dt/ + \X 9x 9y r 8<p y x V Hy *7 


pa 

r 


^av - w 2 cos /3y + vw cos /3^j sin 6 


+ (au - w 2 cos /3x + uw cos /3 Jj cos 6 


(C6) 


where (d/dt) + denotes the total time derivative along the bi characteristic and is written 
as 


I(u ♦ a cos * ( 


— ) = — + — |u + a cos 

W + <“ 


9 1/ \ 9 

V + a cos /3y)— + -Jw + a cos /3 J— (C7) 


The bicharacteristic has the slopes 




For the purposes of the paper, it is desired to express the compatibility relation in 
terms of the independent variables t, X, Y, and $ which are given by equations (7). 
This is accomplished with the aid of expressions (8) which relate the partial derivatives 
with respect to t, x, y, and cp to those with respect to r, X, Y, and >h. First 
consider the first group of terms on the right side of equation (C6). These terms can be 
transformed as follows: 
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X 9x 9y r Bcp X 3X X6 \ b / 8Y 6 0Y r 9$ r 6 V ■* 


9U<p 

9Y 


1 9 Ux G/ o ®Ux 0 ®U(/A I 

+ t( cos ftc — + cos % -5^ + cos Z 3 ^ T^r) + r t~ 


x ax 6 


8Y 


8Y 


r 9$ 


ii^ + S-L(e s .u) + I^ 

x ax 6 9Y ' 1 r a$ 


1 ^^x + 1 

x ax r a$ 


The last step follows because 
e s • U = 0 

Similarly, the expression (C7) for the total derivative with respect to time along the 
bicharacteristic can be transformed as follows: 


— ^ = — — - — + -(u + a cos / 3 x') 
dt/ + 3r 6 8Y x[ ^7 


— — -(s’ + y^)— + — (v + a cos fiy)— 
[ax 6 \ b /9YJ 6\ v 9Y 


+ — w + a cos 
r ' 




= — + -fu + a cos SjA— + - G(V + a) - 6 — + — (w + a cos /3 J) — 
0T A\ /9X 6L J8Y rV 7 9$ 


Thus, the compatibility relation can be written in terms of the independent variables t, 
X, Y, and $ as 


+ — (u + a cos ^x)— + — 
ar XV JdX 6 


G(V + a) - 6 


+I(W + 

8Y r \ 


a cos + Pa<(-7^ + 7(11 + a cos ^x)- 9V 


9$ 
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+ — 
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G(V + a) - + |(w + a cos /3^ = -pa 
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1 8 U X [ 1 8 U^\ | K 
X 9X + * 9$ / + X 
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+ — (av - u^cos /3y 


+ uv cos 


Px) + S1 f e ( av ■ w 2 cos /3y + vw cos /3^ 


+ C °r ^ (a u ~ w^cos /3x + uw cos p^j 


52 



REFERENCES 


1. Rusanov, V. V.: Three-Dimensional Flow About an Arbitrary Blunt Body. Aerospace 

Proceedings 1966, Vol. 1, Joan Bradbrooke, Joan Bruce, and Robert R. Dexter, 
eds., Macmillan Co., 1966, pp. 291-301. 

2. Bohachevsky, Ihor O.; and Mates, Robert E.: A Direct Method for Calculation of the 

Flow About an Axisymmetric Blunt Body at Angle of Attack. AIAA J., vol. 4, no. 5, 
May 1966, pp. 776-782. 

3. Moretti, Gino; and Bleich, Gary: Three-Dimensional Flow Around Blunt Bodies. 

AIAA J., vol. 5, no. 9, Sept. 1967, pp. 1557-1562. 

4. Xerikos, J.; and Anderson, W. A. : A Time-Dependent Approach to the Numerical 

Solution of the Flow Field About an Axisymmetric Vehicle at Angle of Attack. 
DAC-62279 (Contract NAS 8-21141), Douglas Aircraft Co., June 1968. (Available 
as NASA CR-61982.) 

5. Cohen, Gerald A.; Foster, Richard M.; and Dowty, James R.: Synthesis of Optimum 

Structural Designs for Conical and Tension Shell Mars Entry Capsules. NASA 
CR-1365, 1969. 

6. Masson, B. S. ; Taylor, T. D.; and Foster, R. M.: Application of Godunov's Method to 

Blunt-Body Calculations. AIAA J., vol. 7, no. 4, Apr. 1969, pp. 694-698. 

7. Barnwell, R. W.: Three-Dimensional Flow Around Blunt Bodies With Sharp Shoulders. 

AIAA Paper No. 71-56, Jan. 1971. 

8. Barnwell, Richard W.: Time-Dependent Numerical Method for Treating Complicated 

Blunt-Body Flow Fields. Analytical Methods in Aircraft Aerodynamics, NASA 
SP-228, 1970, pp. 177-195. 

9. Barnwell, Richard W.: A Time-Dependent Method for Calculating Supersonic Blunt- 

Body Flow Fields With Sharp Corners and Embedded Shock Waves. NASA 
TND-6031, 1970. 

10. Moretti, Gino; and Abbett, Michael: A Time-Dependent Computational Method for 

Blunt Body Flows. AIAA J. , vol. 4, no. 12, Dec. 1966, pp. 2136-2141. 

11. Brailovskaya, I. Yu.: A Difference Scheme for Numerical Solution of the Two- 

Dimensional, Nonstationary Navier-Stokes Equations for a Compressible Gas. 

Sov. Phys. — Doklady, vol. 10, no. 2, Aug. 1965, pp. 107-110. 

12. Richtmyer, R. D.; and Morton, K. W.: Stability Studies With Difference Equations. 

(I) Non-Linear Stability. (II) Coupled Sound and Heat Flow. NYO- 1480-5 (Con- 
tract AT(30-1)-1480), Courant Inst. Math. Sci., New York Univ., Aug. 25, 1964. 


53 



13. Barnwell, Richard W.: Inviscid Radiating Shock Layers About Spheres Traveling at 

Hyperbolic Speeds in Air. NASA TR R-311, 1969. 

14. Guderley, K. Gottfried: Singularities at the Sonic Velocity. Tech. Rep. 

No. F-TR-1171-ND, Air Materiel Command, U.S. Air Force, June 1948. 

15. Friedman, Manfred P.: Two-Dimensional and Axisymmetric Rotational Flows Past 

a Transonic Corner. J. Aerosp. Sci., vol. 29, no. 4, Apr. 1962, pp. 503-504. 

16. Shifrin, E. G.: On the Extremal Properties of Entropy on a Critical Streamline. 

J. Appl. Math. Mech., vol. 31, no. 3, 1967, pp. 618-620. 

17. Webb, H. G., Jr.; Dresser, H. S.; Adler, B. K.; and Waiter, S. A.: Inverse Solution of 

Blunt-Body Flowfields at Large Angle of Attack. AIAA J., vol. 5, no. 6, June 1967, 
pp. 1079-1085. 

18. Anon.: U.S. Standard Atmosphere, 1962. NASA, U.S. Air Force, and U.S. Weather 

Bur., Dec. 1962. 

19. Kendall, James M., Jr.: Experiments on Supersonic Blunt-Body Flows. Progr. Rep. 

No. 20-372 (Contract No. DA-04-495-Ord 18), Jet Propulsion Lab., California Inst. 
Technol., Feb. 27, 1959. 

20. Stallings, Robert L. , Jr.; and Howell, Dorothy T.: Experimental Pressure Distribu- 

tions for a Family of Blunt Bodies at Mach Numbers From 2.49 to 4.63 and Angles 
of Attack From 0° to 15°. NASA TN D-5392, 1969. 

21. Stallings, Robert L., Jr.; and Tudor, Dorothy H.: Experimental Pressure Distribu- 

tions on a 120° Cone at Mach Numbers From 2.96 to 4.63 and Angles of Attack From 
0° to 20°. NASA TN D-5054, 1969. 

22. Allison, Dennis O.: Calculation of Thermodynamic Properties of Gas Mixtures at 

High Temperatures. M.S. Thesis, Virginia Polytech. Inst., 1965. 

23. Browne, W. G. : Thermodynamic Properties of the Earth's Atmosphere. Radiation 

and Space Phys. Tech. Mem. No. 2, Missile Space Div., Gen. Elec. Co., Nov. 15, 

1962. 

24. Richtmyer, Robert D.: A Survey of Difference Methods for Non-Steady Fluid Dynamics. 

NCAR Tech. Notes 63-2, Nat. Center Atmos. Res., Aug. 27, 1962. 

25. Van Leer, B.: Stabilization of Difference Schemes for the Equations of Inviscid Com- 

pressible Flow by Artificial Diffusion. J. Comput. Phys., vol. 3, no. 4, Apr. 1969. 
pp. 473-485. 

26. Von Mises, Richard: Mathematical Theory of Compressible Fluid Flow. Academic 

Press, Inc., 1958, pp. 100-116. 


54 












(a) Bow-shock-wave and sonic-line locations in plane of symmetry. 

Figure 7.- Comparison of results for equilibrium air flow field about Apollo command 
module at 22° angle of attack and traveling at a speed of 6935 m/s (22 754 ft/sec) 
at an altitude of 45.866 km (150 480 ft). 
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(b) Surface pressure distribution for plane of symmetry (<p = 0° and 180°) and 
plane normal to plane of symmetry ( cp = 90°). 


Figure 7.- Concluded. 
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[oj 
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Distribution line 


O Modified finite-difference method 
□ Unmodified finite-difference method 
A Method of characteristics (ref. 9) 
Experiment (ref. 19) 


(a) Surface pressure distribution. 


(b) Pressure distribution across shock 
at shoulder along line shown. 


Figure 9.- Pressure distribution on surface and along line across shock layer at shoulder for flat-face cylinder 

at 0° angle of attack with Moo = 2.81 and y = 1.4. 










(a) Surface pressure distribution in plane of symmetry. 

Figure 11.- Flow field about spherically blunted cone with sharp shoulder, semiapex 
angle of 60°, and ratio of nose radius to base radius of 0.25 for angles of attack 
of 10° and 20°. = 4.63; y = 1.4. 
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c p , deg 

(b) Peripheral distribution of cross-flow component 
of velocity w for a = 10°. 



<P , deg 

(c) Peripheral distribution of cross-flow component 
of velocity w for a = 20°. 

Figure 11.- Continued. 
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Figure 11.- Continued. 





(f) Surface pressure distribution in plane cp = 90°. 
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(g) Surface distribution of u-component of velocity in plane cp = 90°. 
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(h) Surface distribution of cross -flow component of velocity w in plane cp = 90° 

Figure 11.- Concluded. 
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